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ABSTRACT 

Observational and theoretical evidence suggests that coronal heating is impul- 
sive and occurs on very small cross-field spatial scales. A single coronal loop could 
contain a hundred or more individual strands that are heated quasi-independent ly 
by nanoflares. It is therefore an enormous undertaking to model an entire ac- 
tive region or the global corona. Three-dimensional MHD codes have inadequate 
spatial resolution, and ID hydro codes are too slow to simulate the many thou- 
sands of elemental strands that must be treated in a reasonable representation. 
Fortunately, thermal conduction and flows tend to smooth out plasma gradi- 
ents along the magnetic field, so "OD models" are an acceptable alternative. We 
have developed a highly efficient model called Enthalpy-Based Thermal Evolu- 
tion of Loops (EBTEL) that accurately describes the evolution of the average 
temperature, pressure, and density along a coronal strand. It improves signifi- 
cantly upon earlier models of this type — in accuracy, flexibility, and capability. 
It treats both slowly varying and highly impulsive coronal heating; it provides 
the time-dependent differential emission measure distribution, DEM(T), at the 
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transition region footpoints; and there are options for heat flux saturation and 
nonthermal electron beam heating. EBTEL gives excellent agreement with far 
more sophisticated ID hydro simulations despite using four orders of magnitude 
less computing time. It promises to be a powerful new tool for solar and stellar 
studies. 

Subject headings: hydrodynamics — methods: numerical — Sun: corona — Sun: 
transition region — stars: coronae 



Introduction 



An abundance of observational and theoretical evidence indicates that much of the 
corona is highly dynamic and evolves in response to heating that is strongly time-dependent. 
The evidence further suggests that the cross-field spatial scale of the heating is very small, 
so that unresolved structure is ubiquitous. In particular, many if not all coronal loops are 
bundles of thin strands that are heated impulsively and quasi-randomly by nanoflares. It is 
estimated that a single loop contains several tens to several hundreds of such strands. See 
Klimchukl ( 120061 ) for a detailed justification of these ideas and references to relevant work. 



Three-dimensional (3D) magnetohydrodynamic simulations are extremely useful for 
studying the source of coronal heating (instabilities of electric current sheets, reconnec- 
tion, turbulence, etc.), but they cannot adequately address the complexity that is present 
in a single coronal loop, much less an entire active region. A more feasible approach is to 
treat the magnetic field as static and to solve the one-dimensional (ID) hydrodynamic equa- 
tions along many representative flux strands using an assumed heating rate. The individual 
strands must be treated separately. It is not valid to approximate a loop as a monolithic 
structure with uniform heating corresponding to the average for the component strands. 
This gives a completely different and incorrect result. 

There is reason to believe that the diffuse cor ona that lies between d istinct bright 
loops is also comprised of elemental strands (e.g., lAschwanden et al.l 120071 ). If roughly 



100 strands are present in a single loop, then the numbers present in active regions and 
the global Sun are truly staggering. While it is possible t o construct a limited numbe r 
of model active regions with time-dependent ID simulations (IWarren &: Winebargerl 120071 ). 
it is not possible to investigate a wide range of values for the coronal heating parameters 
that must be assumed, such as the depend ence on magnetic field strength, loop length, 
etc. ( iMandrini. Demoulin. fc Klimchukl |2000| ). This is a major limitation, since we are still 
struggling to identify the properties and physical origin of the heating mechanism. Progress 
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in the foreseeable future must therefore rely on simplified solutions to the hydro equations 
that treat field-aligned averages and are much less computationally intensive. These are 
sometimes called "OD models" because there is only one value of temperature, pressure, and 
density at any given time in the simulation. 



topics, including c oronal l oops (Cargil l fc K 



OD models have been developed p reviously bvlFisher fc Hawleyl (119901 ) and lKopp fc Poletto 
( 1993 ). but the best known is t hat of Cargilll (1994 ). It has been used to study a variety of 
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2006 



Klimchuk & Cargilll 2001; 



Parenti et ap2006h . fl ares JReeves fc WarreJ2002l ; |Patsourakos. Antiochos. fc Khmchukll2002ri 



post- eruption arcades (jReeves fc Forb esl l2005h . and active stellar coronae JCargill fc Klimchukl 
20061 ). We have learned a great deal with the Cargill model, and our understanding has 
now advanced to the point where a more accurate and flexible model is required. This ar- 
ticle presents an improved 0D model called Enthalpy-Based Thermal Evolution of Loops 
(EBTEL). As the name suggests, a key aspect of the model is an explicit recognition of the 
important role that enthalpy plays in the energy budget. 

EBTEL improves upon the Cargill model in several important ways. First, whereas the 
Cargill model is limited to an instantaneous heat pulse, EBTEL accommodates any time- 
dependent heating profile and can include a low-level background heating if desired. Second, 
EBTEL accounts for thermal conduction cooling and radiation cooling at all times during 
the evolution. The Cargill model assumes that only one or the other operates at any given 
time. Third, EBTEL has options for heat flux saturation and nonthermal electron beam 
heating. Finally, EBTEL is unique among 0D models in that it provides the time-dependent 
differential emission measure distribution of the transition region footpoints. Emission from 
the transition region plays a critical role in spatially unresolved observations, such as stellar 
observati ons and ob servations of the solar spectral irradiance, which is important for space 
weather (jLeanl 119971 ). Note that footpoint emission is not limited to the cooler (< 1 MK) 
plasma traditionally associated with the transition region. It can also include hot emissions 
that originate from the base of very hot loops. The so-called moss seen in the "c oronal" chan- 



nels of the Transition Region and Coronal Explorer (TRACE) is an example (IBerger et al. 
19991 : iMartens. Kankelborg. fc Bergerlboooh . 



We describe the coronal and transition region parts of EBTEL in the next two sections. 
We then present example simulations and compare them with corresponding simulations 
from a ID model and, in one case, the Cargill model. We conclude with a discussion of 
EBTEL and the possible significance of the example simulations. 
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2. Coronal Model 

Both OD and ID models of the corona are traditionally referred to as "loop" models. 
However, as we have discussed, there is good reason to believe that the loop-like intensity 
features seen in coronal images are actually comprised of many individual strands. EBTEL 
treats these individual strands, which are mini flux tubes in which the plasma is approxi- 
mately uniform within the cross section. Multiple strand models can be combined to form a 
loop bundle. 

Under static equilibrium conditions, the coronal portion of a strand is characterized 
by an exact balance bet ween energy input (coronal heat i ng) and energy losses by radiation 



and thermal conduction (jRosner. Tucker ] &: Vaianalll978l ; ICraig. McClymont. &: Underwood 
1978 ; Vesecky. Antiochos. &: Underwood 19791 ). Some of the coronal heating energy — less 



than half — is radiated directly to space, and a heat flux carries the remainder to the transition 
region, from where it is more efficiently radiated. Temporal variations in the heating rate 
produce a well-defined response involving the transfer of mass between the chromosphere 
and corona. Heating increases cause the coronal temperature to rise and the downward heat 
flux to intensify. The transition region is unable to radiate the extra energy, and heated 
plasma flows into the strand in response to enhanced pressure gradients. This is the well- 
know process of chromospheric "evaporation." The inverse process ("condensation") occurs 
when the heating rate decreases. As the coronal temperature declines, the reduced heat flux 
is insufficient to power the transition region radiation. The plasma cools, pressure gradients 
drop to sub- hydrostatic values, and material drains from the strand. 

The basic idea behind EBTEL is to equate an enthalpy flux of evaporating or condensing 
plasma with any excess or deficit in the heat flux relative the transition region radiation 
loss rate. An excess heat flux drives evaporative upflows, while a deficient heat flux is 
compensated by condensation downflows. The key assumption of the model is that the 
radiative losses from the transition region and corona maintain a fixed ratio at all times. This 
ratio is the same one that applies during static equilibrium conditions. We defer justification 
of this assumption until later and now derive the equations that define the model. 

We begin with the ID time-dependent equation for energy conservation: 

arp a a an 

where 

E = -P+-pv 2 (2) 

is the combined thermal and kinetic energy density; s is the spatial coordinate along the 
magnetic field; n, T, P, and v are the electron number density, temperature, total pressure, 
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and bulk velocity, respectively; F is the heat flux; Q is the volumetric heating rate; g» is 
the component of gravity along the magnetic field; and A(T) is the optically thin radiative 
loss function, for which we use a piecewise continuous form based on the atomic physics 
calculations of J. Raymond (1994, private communication) and twice the coronal elemental 
abundances of Meyer ( 19851 ): 
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(3) 



The highest tem perature range of the loss function is dominated by thermal bremsstrahlung 
(jCox fc Tuckerlll969l ). Equation ([I]) assumes a constant cros s-sectional area, which is appro- 
priate for distinct coronal loops and their constituen t strands (lKlimchukll2000l ; IWatko and Klimchuk 
20001 ; lLopez Fuentes. Klimchuk. fc Demoulinll2006l ). but probably not for the diffuse corona. 
We also assume that the loop is symmetric, so only one half need be considered. We define s 
to increase from footpoint to apex. The downward heat flux is therefore a negative quantity. 
To simplify the discussion, we do not at this time include the energy and particle fluxes of 
a possible nonthermal electron beam. These will be added later. 



If the flow is subsonic (v < C s = 1.5 x lO^ 1 / 2 = 2.6 x 10 7 cm s" 1 at T = 3 MK) and 
the loop is shorter than a gravitational scale height (z apex < H g = 5.0xl0 3 T = 1.5xl0 10 cm 
at T = 3 MK), then the kinetic energy and gravity terms in equation ([1]) can be neglected, 
leaving 



5 d 



(Pv) 



OF 



Q - n 2 A(T). 



(4) 



3dP 

2~dt 2ds v ~ ds 
We now define the base of the corona, designated by subscript "0," to be the location where 
thermal conduction changes from being a cooling term above to a heating term below. This 
occurs at the top of a thin transition region, close to the chromospheric footpoint. Integrating 
equation (T4j) over the coronal portion of the strand and noting that the velocity and heat 
flux both vanish at the apex due to symmetry, we obtain 



3 dP 5 

-L— w -P v + F + LQ — TZ C 
2 at 2 



(5) 



The over bars indicate spatial averages along the coronal section, which has length L from 
the coronal base to apex. The first two terms on the right-hand-side of equation (jSJ) are 
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the enthalpy flux and heat flux at the coronal base, and 7Z C is the radiative cooling rate 
per unit cross sectional area in the corona (erg cm -2 s _1 ). Since temperature, pressure, and 
density typically vary by less than a factor of 2 along the coronal section, the averages are 
quite characteristic of the entire section. We sometimes refer to them as simply the coronal 
values. 

If we instead integrate equation (j4j) over the transition region, spanning from the top of 
the chromosphere to the base of the corona, we obtain a similar result: 

2 l ~&t ~ ~2 P ° V ° ~ Fo + lQtr ~ Utr ' (6) 

except that the spatial averages are now along the transition region, which has thickness / 
and radiative cooling rate lZt r . In deriving this result, we made use of the fact the heat flux 
and enthalpy flux are both ignorable at the top of the chromosphere. During evaporation, 
a very small heat flux does in fact reach the top of the chromosphere, but most of the heat 
flux is dissipated throughout the transition region, heating each layer to the next higher 
temperature. No heat flux reaches the chromosphere during condensation. 

Concerning the enthalpy flux at the top of the chromosphere, we note that conservation 
of mass requires that the electron flux be nearly constant through the transition region 
during both evaporation and condensation: 

J = nv « J . (7) 

Together with the ideal gas law, 

P = 2knT, (8) 

where k is Boltzmann's constant and we have assumed a fully ionized hydrogen plasma, 
equation (J7]) implies that the enthalpy flux is proportional to temperature. The enthalpy 
flux is therefore much smaller at the top of the chromosphere than at the base of the corona 
and can be safely ignored in equation (jSJ). 

Because the transition region is so thin, we can neglect the terms involving I in equation 
([6]) and are left with 

^P v « -F - TZ tr . (9) 

When |i*o| > TZtr, there is an excess heat flux that drives a positive enthalpy flux (evapora- 
tion). When |Fo| < 7lt r , there is a negative enthalpy flux (condensation) that combines with 
the heat flux to power the radiation. Static equilibrium corresponds to an exact balance 
l-^o I — TZ tr . 
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Combining equations (jSJ) and ([9]), we obtain the following equation for the evolution of 
the coronal pressure: 

dP 2 



dt 3 



Q-\in c + n tr ) 



(10) 



We note that the same equation is obtained if we keep the terms involving I in equation ([6]) 
and interpret P and Q as the spatial averages along the entire strand, including both the 
transition region and coronal sections, with L then being the total length. Equation (ITUl) 
reflects the energetics of the combined corona-transition region system. Energy enters the 
system only through coronal heating, and energy leaves the system only through radiation. 
Thermal conduction and flows transport energy between the corona and transition region, 
but they do not add or remove energy from the system. 

It has been suggested that most "coronal" heating occurs in the transition region, in 
which case lQ tr >> LQ and equation (j5J) is not a good approximation. We do not believe 
this is a likely possibility, however. For one thing, the transition region is very thin. For 
another, it moves up and down a significa nt distance in re sponse to changes in the spatially 



integrated heating rate within the strand (jKlimchuk!l2006l ) . One might expect the positional 



dependence of the heating to be more closely related to the magnetic field than to the variable 
location of the transition region plasma. It is nonetheless possible to use EBTEL to study 
direct heating of the transition region. Its effect is very similar to that of a nonthermal 
electron beam, which is discussed in Sections 2.1 and 4.5. A minor difference is that an 
electron beam will slightly decrease the coronal mass. 

We have defined the transition region to be the section of the strand where the heat flux 
is an energy source term. By this definition, its thickness is roughly 10% of the strand half 
length (l/L ~ 0.1). This is not exceptionally small. However, most of the coronal heat flux is 
deposited within the extreme lower part of transition region, which is also where most of the 
radiation is emitted. We could therefore redefine the transition region to be much thinner 
and our model would be substantially unchanged. As an example, consider an equilibrium 
strand with an apex temperature of 2 MK and half length L = 7.5 xlO 9 cm. 80% of the heat 
flux is deposited over a distance of only l/L = 0.013, and 50% of the heat flux is deposited 
over an even shorter distance of l/L = 0.00069. 

We wish to express the basic pressure equation, equation (jTUJ), in terms of the time- 
dependent variables P, T, and n. We therefore approximate the radiative loss rate from the 
corona as 

TZ c ^n 2 A(f)L. (11) 

This would be exact if the coronal density and temperature were perfectly uniform instead 
of approximately so. 
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Next, we assume that the radiative loss rates of the transition region and corona maintain 
a fixed ratio at all times: 



Since we want the model to apply during slow evolution as well as fast, C\ should be equal to 
the static equilibrium value. One difficulty is that 1Z tr /1Z c is different for different equilibrium 
conditions. In particular, it depends on the apex temperature of the strand, T a . Table [1] 
lists lZt r /lZ c determined from exact equilibrium solutions in a semi-circular strand of half 
length L = 2.5 x 10 9 cm. Six apex temperatures ranging from 0.8 to f 0.4 MK correspond 
to six different spatially- uniform heating rates. Except for lowest temperature case, TZtr/T^-c 
increases monotonically with T a from 1.8 to 20.7 MK. In one implementation of EBTEL, 
we let c\ vary according to a third order polynomial fit to these data. However, after some 
experimentation, we found that a constant value c\ = 4.0 provides the best overall agreement 
with ID simulations, especially in cases of impulsive heating. Table [2] lists Ht r /1Z c for a 
longer equilibrium strand with L = 7.5x 10 9 cm. The ratio is reasonably close to 4 for apex 
temperatures ranging from 1-4 MK. All of the results presented in this paper use a constant 
value c\ = 4.0. 

We do not yet have a compelling physical argument for why TZ tr /H c should be constant 
even when the strand is far from equilibrium. Fortunately, this does not appear to be an 
important assumption, at least not for the simulations presented in this paper. It is certainly 
not important during times of strong evaporation, when the evolution is essentially a balance 
between the downward heat flux and upward enthalpy flux, and radiation plays no significant 
role. The radiative losses are only 10 -3 of the heat flux during the strong evaporation phase 
of the nanoflare simulation of Section 4.1 (example 1). Radiation is very important during 
times of strong condensation, on the other hand. The assumption that 1Z^ r /TZ c equals the 
equilibrium value could then potentially cause problems. It turns out that condensation 
is fairly mild in all of our example simulations in the sense that the radiative losses never 
greatly exceed the heat flux. They are at most a factor of 3.6 larger in example 1. 

We have considered a wide variety of heating scenarios, some discussed in Section 4 
and others not reported in this paper. The fact that EBTEL is able to reproduce exact ID 
solutions as well as it does gives us considerable confidence in the approximations of the 
model, including the assumption that lZ tr /lZ c = 4.0. 

With equations ([IT]) and f|T2|) . we can now express equation ffTUj) for the evolution of 
coronal pressure in terms of the fundamental variables P, T, and n. 

We next move on to an equation for the coronal density. The total mass contained in 
the coronal section of the strand changes as material evaporates and condenses. Specifically, 
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the time derivative of the electron column density nL (electrons per unit cross sectional area) 
is equal to the flux of electrons through the coronal base: 

§- f (nL) = J . (13) 

This can be derived trivially by integrating the ID equation of mass conservation from the 
base of the corona to the apex. Combining equations ((71), (jHJ), and Qj, we can write the 
electron flux as 

Jo = ~5^ (F ° + 7 ^ r) - (14) 
Substituting into equation f[T3"l) . we get 

dn Co , „ m x / % 

A=-5^ (F ° + K "»' (15) 

where we have introduced C2 for the ratio between the average coronal temperature and apex 
temperature, 

C2 = (16) 
and C3 for the ratio between coronal base temperature and apex temperature, 

ca = (17) 

Tables [1] and [2] list the values of these two ratios for the exact equilibrium solutions in the 
short and long strands, respectively. T/T a is very close to 0.87 in all cases, while T /T a varies 
over a fairly narrow range from 0.22 to 0.61. In the implementation of EBTEL with variable 
Ci, we also let C3 vary based on a polynomial fit to the data in Table [TJ As already indicated, 
constant values give better overall agreement with the ID simulations; C3 = 0.5 seems to 
work best and is the value used in the examples presented here. Note that C2 = 0.87 is not 
far from 7/9, which corresponds to a constant heat flux solution. 

Using equations (11 II) and (TT2"j) . we can express TZ tr in terms of our fundamental variables 
n and T, but we still need an expression for Fq. The classical expression for the heat flux is 

F c = -*oT 5 / 2 ^, (18) 
where k = 1.0 x 10~ 6 in cgs units. Noting that 

r5/2 f - ?! ^ • (19) 

we can approximate the heat flux at the base as 

Fc«~Ko-^-, (20) 
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where T a = T/c?. The precise value of the coefficient depends on the details of the temper- 
ature profile; 2/7 corresponds to a constant heat flux, while 4/7 corresponds to a constant 
heat flux divergence. 

The classical heat flux is, however, unphysically large during times of exceptionally 
high temperature and/or exceptionally low density, such as during the earliest phase of an 
impulsive heating event. Under these conditions, the heat flux saturates at approximately 

Fs - -4^ T3/2 > (21) 



wher e m e is the electron mass a nd (3 is a flux limiter constant that we set to 1/6 (jLuciani. Mora, fc Virmont 



19831 ; iKarpen fc DeVordll987l ). We consider two possibilities in our simulations. First, we 

brm 

(22) 



set F = F c at all times, regardless of the temperature and density. Second, we use the form 

F r F« 



which reduces to F c when \F C \ -C \F S \ and to F s when \F C \ 3> \F S \. We can now express 
equation ( |T5"i) for the coronal density evolution in terms of the fundamental variables. 

The last governing equation, for the coronal temperature evolution, follows straightfor- 
wardly from the ideal gas law: 

df nfldP ldn\ . , 



dt \P dt n dt J 

Note that this is not exact because the ideal gas law is not exact when average values of P, 
T, and n are used. 

Summarizing, the coronal part of EBTEL is defined by: evolutionary equations fllOp . 
( TT5|) . and fl23l ; the assumption given by equation (fT2l) ; the approximations given by equations 
f fTT]) and fl20l) or fl22l) ; and parameters c\ = 4.0, C2 = 0.87, and C3 = 0.5. Note that C2 and C3 
appear only together as a ratio in equation f|T5|) . so there are really only two parameters in 
the model. 

The plasma velocity at the base the corona can be obtained straightforwardly from the 
electron flux, equation f|T4l) . according to 

c 3 2kTJ 

vo = - — . (24) 
Using T in place of T gives the velocity at that temperature in the transition region. 
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2.1. Nonthermal Electron Beam 



The mechanism that directly heats the coronal plasma may also produce energetic par- 
ticles. It is thought, fo r example, that a sizable fraction of the total en ergy of a flare goes into 
nonthermal electrons (jSaint-Hilaire fc Benzll2002l ; lEmslie et al.ll2005l ). We therefore have in- 
corporated a nonthermal electron beam into EBTEL. We assume that the electrons originate 
from the existing strand plasma and stream freely along the magnetic field to the coronal 
base. We further assume that all of the beam energy goes into the enthalpy of evaporating 
plasma and that any chromospheric radiation that may be produced is negligible. Because 
we do not consider the details of the energy deposition (i.e., how it depends on column 
depth), our calculation of the differential emission measure of the transition region (Section 
3) is not reliable when nonthermal electrons are included. 

The effect of the electron beam on the coronal energy budget is straightforward. The 
corona gains energy from the enthalpy of the evaporated plasma, but it loses energy because 
electrons must be removed from the thermal pool to supply the seed particles for the beam. 
In general, the gain far exceeds the loss because the mean energy of the accelerated electrons, 
£, is much greater than their thermal energy, (3/2) kT. 

If T and J are the energy flux and particle flux of the beam, respectively, so that 



T = £J, 



(25) 



then we must modify our equations by subtracting T from the right of equation 0, adding 
{3/2)kTJ to the right side of equation §5§, and adding J to the right side of (TT31) . Note 
that T and J are both negative quantities. The evolutionary equations for pressure and 
density are then 

%\q 



dP 
~dt 



-jr [He + Ktr) ~ j 



and 



dn c 2 JF 

Tt = -^kUT [Fo + ntr) + TL 



1 



3 kT 
Y~£ 

c 2 £ 
5c^kf 



(26) 



(27) 



We have avoided the tricky issue of electron return currents. Quasi neutrality of the 
plasma requires either that protons are accelerated with the electrons, which is thought to 
be unlikely, or that an electron return current replenishes the electrons lost to the beam. 
With a return current, the unity term inside the last set of parentheses disappears from 
equation (j2"Tj) . and the temperature in the kT JS term of equation fl26|) must be replaced by 
the temperature difference between the strand plasma and the replenishing electrons. Note 
that both of these terms are negligible as long as £ 3> kT, even without a return current. 
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2.2. 



Differential Emission Measure 



Most observed plasmas are multi thermal, even within a single observational pixel. An 
important quantity is therefore the differential emission measure, DEM(T), which describes 
how the plasma is distributed in temperature. Spatial variations in the coronal temperature 
tend to be greatest across the magnetic field. In a multi-stranded loop bundle, for example, 
the different strands will have different temperatures if the heating is steady but unequal or 
if it is impulsive but out of phase. There is also some temperature variation along the field. 
We here consider the differential emission measure of a single stand of unit cross sectional 
area: 



The transition region is treated carefully in the next section. For the corona, we make the 
crude approximation that the total emission measure, 2Ln 2 , is distributed uniformly over 
the temperature interval 0.74T a < T < T a . The average temperature T falls exactly in the 
middle of this interval. Note that our approximation is not critical, because the DEM(T) of 
a coronal observation is determined primarily by the distribution of different strands rather 
than by the variation along each strand. 



The situation is very different in the transition region, where the temperature and 
density vary dramatically over a short distance along the magnetic field. Emission from the 
transition region is a critically important component in many observations, such as spatially 
unresolved observations of stars or measurements of the full Sun spectral irradiance. Even 
high-resolution observations on the solar disk tend to have lines of sight that pass through 
both coronal and transition region plasmas. We do not attempt to model the detailed spatial 
structure of the transition region, but instead deal directly with the differential emission 
measure. We have developed two separate approaches. The one we now discuss is the 
easiest to implement and has been used for all of the examples shown in the paper. The 
second approach, presented in the Appendix, has the advantage of being physically more 
revealing. It treats the limiting cases of strong evaporation, strong condensation, and static 
equilibrium, and provides simple analytical expressions for DEM(T) in each case. The two 
approaches produce similar results. Neither is valid when nonthermal particles are important. 

We begin with the steady state version of the energy equation, equation (jlj), in the 




(28) 



3. Transition Region Model 
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absence of local heating: 

k pi ftp 

2S (P " )+ a7 + n2A(r)=3 °- (29) 

which should be approximately correct in the transition region. We next assume that the 
heat flux term can be approximated by —K T 3 ^ 2 (dT/ds) 2 . This is strictly valid when the 
scale lengths of temperature and heat flux are the same. Rewriting the enthalpy term using 
the ideal gas law and constant mass flux, and assuming that pressure is the same in the 
corona and transition region, the energy equation becomes 

^ 3/2 (f) 2 - 5 *4H^y A(r) ~ o ' <3o) 

This is quadratic in dT/ds and can be solved trivially. DEM(T) then follows directly from 
equation (1251) . 

We have computed the radiative loss rate from the transition region by integrating the 
product DEM(T) x A(T) over the transition region temperature interval and find that it is 
similar to the value obtained from 1Z tr = C\TZ C . The only exception is during times of strong 
evaporation, when the C\ assumption is unimportant. 



Results 



4.1. Example 1 



We have coded up EBTEL in the Interactive Data Language (IDL) and now examine 
several simulations that were run on a desktop computer. The first example considers an 
impulsive energy release in a static equilibrium strand of half length L = 7.5 x 10 9 cm. An 
average coronal temperature of 0.52 MK in the initial equilibrium is produced by a heating 
rate of 10~ 6 erg cm -3 s" 1 . We obtain the equilibrium by guessing at the values of T, n, and P 



using scaling law theory (IRosner. Tucker, fc Vaianalll978l ; ICraig. McClymont. &: Underwood 
19781 ) and letting the strand evolve while holding the heating rate constant. T changes very 
little during the relaxation, while n and P decrease by about a factor of two. 

We impose a nanoflare energy release on top of the steady background heating. It has 
a triangular profile with a total duration of 500 s and a peak value of 1.5 x 10 -3 erg cm -3 
s -1 , 1500 times stronger than the background. Nonthermal electron beams are excluded 
from all but the last our examples. The solid curves in Figure [p show how T , n, and P 
respo nd to the event. The generic behavior is well documented ( ICargilll Il994l ; iKlimchuk 



20061 ). Temperature and pressure rise abruptly as the nanoflare energy is converted into 
thermal energy at a roughly constant density. An intense heat flux drives strong evaporation 
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and the strand begins to fill with plasma. The temperature then declines as the nanoflare 
shuts off, but evaporation continues, and the peak density is not reached until well after the 
nanoflare has ended. Radiation becomes progressively more important as the temperature 
falls and density rises. It eventually takes over from thermal conduction as the dominant 
cooling mechanism. The strand then enters a long phase of draining and condensation. 

We have run an exactly corresponding simulation with our sophisticate d ID hydro code 



called the Adaptively Refined Godunov Solver (ARGOS). As described in lAntiochos et al. 



(119991 ). the code uses an evolving numerical mesh to resolve steep gradients wherever they 
may occur. We use the same radiative loss function used in EBTEL. For the ID simulation, 
we make the additional assumptions, not required of EBTEL, that the strand is semi-circular, 
lies in a vertical plane, and is heated in a spatially uniform manner. We earlier defined the 
boundary between the corona and transition region to be the location where the divergence 
of the heat flux changes sign. This is not practical in the ID simulation due to the more 
complicated temperature structure associated with waves and even shocks that are excited 
by the impulsive energy release. We therefore compute coronal averages by averaging over 
the upper 80% of the strand. These averages are indicated by dashed lines in Figure [U The 
small wiggles are due to the aforementioned waves. 

There is good agreement between the EBTEL and ID results. This is highly encouraging 
given that EBTEL requires approximately four orders of magnitude less computing time. 
This run took only about 10 s. The biggest differences are in the density and pressure, 
where the EBTEL values are about 20% too high for the first 2000 s. The way that the ID 
averages are computed is a contributing factor, since density and pressure are highest in the 
lower part of the strand leg that is excluded from the averages. Another contributing factor 
is that EBTEL assumes that all plasma energy is thermal (f-P)- In fact, some of the energy 
is kinetic, so the pressure is artificially inflated. 

Figure [2] shows the differential emission measure distribution for the full strand averaged 
over the first 10 4 s of the simulation. Both the corona and transition region are included. 
One reason for averaging over time is to simulate the observation of a multi-stranded loop. 
If the strands are heated randomly, then the time average of a single strand is equivalent 
to an instantaneous snapshot of an unresolved bundle. As long as the strands get reheated 
after they cool, then at all times there exists one strand in the bundle for each small time 
interval from the full simulation. As in Figure [TJ EBTEL is represented by the solid curve 
and the ID model is represented by the dashed curve. The agreement is once again very 
encouraging, especially considering that the DEM(T) spans more than 3 orders of magnitude. 
The EBTEL values are too high by factors of 2-3 at the higher temperatures. This is partly 
because of the enhanced densities discussed above and partly because temperature decreases 



- 15 - 



more gradually below 3 MK in the EBTEL simulation (see Figured]). The slower cooling 
rate could be because our assumption of a constant c\ = 4.0 is not quite correct. We are 
currently investigating this issue. 

The DEMiT) plotted here and defined in equation (I2"8j) indicates the amount of plasma, 
n 2 As, that is present in temperature interval AT. It has units of cm -5 K _1 . Some authors 
instead use i 

DEM ln {T) = n 2 (^f) , (31) 

which indicates the amount of plasma present in the logarithmic temperature interval A In T 
and has units of cm -5 . The two definitions differ by a factor T: DEMi n (T) = TxDEM(T). 
Figure [3] shows DEMi n (T). Still other authors define the differential emission measure in 
terms of the base 10 logarithm: DEM log (T) = ln(10) x Tx DEM(T). 

Figure H] separates out the contributions to DEM(T) from the coronal (dashed) and 
transition region (dot-dashed) sections of the EBTEL simulation and the coronal section of 
the ID simulation (dotted). The transition region contribution of course dominates at low 
temperatures, but it is also significant at hotter temperatures that are normally associated 
with the corona. The transition region and coronal contributions are equal at T = 1.0 MK, 
which is approximately 1/4 the temperature of the hottest significant emission measure. 
Although it is difficult to observe the transition region in isolation from the corona, since 
lines of sight that reach the transition region must pass through the corona, it is easy to 
observe the corona in isolation from the transition region simply by looking above the limb. 
The agreement between the coronal DEM(T) curves from the EBTEL and ID simulations 
is reasonably good, except below 0.5 MK, where the EBTEL values are unreliable (since the 
DEM is rather arbitrarily cut off at 0.74T a ; Section 2.2). We discuss the agreement with 
actual observations in Section 5. 

The results presented above assume a classical heat flux at all times. We have repeated 
the simulations with a heat flux that is allowed to saturate according to equation (122]) . 
The results differ only early during the nanoflare energy release, when saturation limits the 
thermal conduction cooling and the average coronal temperature rises to maxima of 4.7 MK 
and 7.6 MK in the EBTEL and ID simulations, respectively (peak apex temperatures are 
of course higher). The densities are very low at this time, however, so the time-averaged 
DEM(T) is minimally affected. The curves are nearly indistinguishable from those in Figure 
[2] the only difference being that the high temperature tail is extended by about 0.02 in the 
logarithm. We note that very hot emission, though generally very faint, provides the best 



diagnostics of nanoflare properties (jPatsourakos fc Klimchuki 120061 ) . Care should be taken 



to include heat flux saturation when studying the hottest emission. 
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4.2. Example 2 

We next consider a much more impulsive nanoflare. It has the same total energy as the 
first example (5.625 x 10 9 erg cm -2 ), but the duration is ten times shorter (50 s) and the 
amplitude is ten times larger. This scenario provides a much better comparison with the 
Cargill model, in which all of the energy is deposited instantaneously. Figure [5] shows the 
evolution of T, n, and P as given by EBTEL (solid), the ID model (dashed), and the Cargill 
model (dotted). Only classical heat flux results are presented, since the Cargill model does 
not include saturation effects. EBTEL again reproduces the ID results quite well, although 
the densities and pressures are about 50% too high over the first half hour. 

The EBTEL and Cargill results differ in several important respects. T reaches a maxi- 
mum of 8.1 MK in the EBTEL model, similar to 7.7 MK in the ID model, whereas the Cargill 
model peaks at only 4.0 MK. The Cargill model predicts substantially higher densities and 
pressures throughout most of the simulation. The primary reason is the assumption that 
radiation is ignorable during the first phase of cooling (ending at 1700 s). Since radiation 
is the only mechanism by which energy can leave the system, the thermal energy density 
and therefore the pressure are constant. Another reason for the excess pressures in both 
the Cargill and EBTEL models is the neglect of kinetic energy. All of the plasma energy is 
assumed to be thermal. This is reasonable only when the Mach number is small. The Mach 
number is mostly less than 0.15 after 500 s in the ID simulation, but there are locations in 
the strand where it approaches 3 shortly after the nanoflare ends. 

A final difference in the Cargill model is the prediction of a catastrophic cooling late 
in the evolution, at approximately 8000 s. This is not present in either the EBTEL or ID 
simulations and is a consequence of the fact that no background heating is possible in the 
Cargill model. The radiative loss function, A(T), is such that a thermal instability causes 
the temperature decline to accelerate in the Cargill model until a pre-set limit is reached 
(usually 0.1 MK). In the EBTEL and ID models, the temperature asymptotically approaches 
the static eq uilibrium value co r respo nding to the background heating rate. We note that the 



0D model of iFisher fc Hawleyl (119901 ) predicts a catastrophic cooling even in the presence of 



background heating, but this appears to be a spurious result, at least in some cases. 



4.3. Example 3 

As a third example, we consider a qualitatively different heating scenario. The strand 
begins in static equilibrium with a uniform heating rate of 2xl0 -4 erg cm -3 s _1 . The heating 
rate is slowly reduced by a factor of 100 over a period of 50,000 s, as shown at the bottom 
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of Figure [6j It is maintained at the reduced level for 5000 s, then suddenly increased to the 
original level over 100 s. It is maintained at that level for 3900 s, then suddenly decreased 
again over 100 s. It remains at the reduced level for the remainder of the simulation. 

The top two panels of the figure show the evolution of temperature and density for 
EBTEL (solid) and the ID model (dashed). The 0D solution tracks the ID solution very 
well. Temperature is systematically high, but the detailed shapes of both the temperature 
and density profiles are faithfully reproduced. This shows that our assumption c\ = 4 is 
reasonable for situations other than impulsive heating. 

4.4. Example 4 

The final two examples are modifications of example 1. The 500 s nanoflare is ten times 
more intense in example 4. Figure [7] shows the time-averaged DEM(T) curves for the whole 
strand (solid), corona (dashed), and transition region (dot-dashed). Note that the coronal 
curve is strongly peaked near 3 MK, as observed in active regions, which we return to shortly. 

4.5. Example 5 

Example 5 differs from example 1 only in the form of the nanoflare energy release. Half 
of the nanoflare energy is assumed to go into direct plasma heating, and the other half is 
assumed to go into nonthermal electrons with a mean energy of 50 keV. Figure M shows 
the coronal DEM(T) curve (solid) together with the corresponding curve from example 1 
(dashed). The curves are nearly identical except that plasma hotter than 3 MK is missing 
from example 5. The reason for this difference is as follows. The amount of evaporated 
material is determined largely by the total energy that is released, regardless of its form. 
Temperature, on the other hand, depends strongly on the form of the energy release. With 
direct plasma heating, the coronal temperature rises until either the nanoflare ends or the 
downward heat flux balances the nanoflare heating rate. In contrast, nonthermal electrons 
have no direct effect on the coronal temperature. Note that direct heating of the transition 
region, discussed after equation QUI) , has a similar effect on the coronal evolution as does a 
nonthermal electron beam. 
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4.6. Additional Tests 



We have tested EBTEL against two other ID hydro codes and found good agree- 
ment in both cases. F. Reale kindly simulated example 1 using the Palermo- Harvard code 



(IPeres et al.lll982l ). and K. Reeves kin dly simulated a loop -top flare with a peak temperature 
of 29 MK using the NRLFTM code iMariska et al.lll982h . It is interesting that the plasma 
evolution is similar even though the Palermo and NRLFTM codes use different radiation 
loss functions than do EBTEL and ARGOS. This shows that the precise form of the loss 
function is not important whenever the heating is impulsive. 



5. Discussion 



As evidenced by these examples, our simple OD model is an excellent proxy for more 
sophisticated and far more co mputationally i n tensive ID hydro simula tions. It improves sub- 
stanti ally on the OD models of lCargilll |l99J ) , [Fisher & Hawlevl (jl990l ). and Kopp & Poletto 
(119931 ) . The Cargill model assumes that heating is instantaneous and that cooling occurs 
either by thermal conduction or by radiation, but not by both at the same time. The Fisher- 
Hawley model: (1) predicts abrupt evolutionary changes as the strand evolves between three 
distinct regimes; (2) does not account for the evaporation that continues well beyond the end 
of an impulsive heating event; and (3) cannot return to the pre-event state due to unphysical 
catastrophic cooling. The Kopp-Poletto model shares some similarities with EBTEL, but 
it treats the flows in a fundamentally different way. Like EBTEL, it equates the enthalpy 
carried by evaporative upflows with an excess heat flux, but the excess is determined relative 
to the pre-event state, rather than to the time-varying radiative losses from the transition 
region. Condensation downflows in the model are given by a density-dependent fraction of 
the free-fall velocity. In actuality, gravity plays no direct role in condensation, since the 
downflows are driven by pressure gradient deficits relative to hydrostatic equilibrium, in the 
same way that evaporative upflows are driven by pressure gradient excesses. Gravity sets the 
value of the hydrostatic gradient, but it is only the deficit or excess relative to this value that 
is important for the flows. Inclined strands experience essentially the same condensation and 
evaporation as do upright strands of the same length. Finally, EBTEL has advantages over 
all three of the other models in that it provides the DEM(T) of the transition region and 
treats nonthermal electron beams and heat flux saturation. 



One obvious application of EBTEL is to investigate the idea that the basic structural 
elements of the corona are very thin, spatially unresolved magnet ic strands t hat a re heated 
impulsively. Loops may be bundles of such strands, as reviewed in iKlimchukl (120061 ). and the 
diffuse corona may be similarly structured. Differential emission measure distributions are 
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one important test of this idea. Observed DEM(T) curves from active regions and the quiet 
Sun tend to be peaked near 10 6,5 and 10 61 K, respectively, and to have a slope (temperature 



power law index) > 0.5 coolward of the peak ([Raymond fc Doyldll98ll ; iDere &: Masonlll993 



Brosius et al.lll996l ). This is consistent with the coronal DEM(T) curves of examples 1 and 4 
(Figures H] and [7]). The full loop curves are discrepant, on the other hand, due to the strong 
contribution from the transition region. The cited observations were made on the disk and 
should in principle include the transition region component. However, it is possible that ab- 



of transition region lines used to construct the DEM(T) curves (e.g., 


Daw. DeLuca. & Golub 


2005 




DePontieu et al. 


1999: 


Doschek & Feldman 


1982 




Schmahl & Orrall 


1979). We are 



currently investigating the magnitude of this effect. 

One of the great mysteries of coronal physics that has come to light in the last few 
years is the discovery that warm (~ 1 MK) coronal loops are much denser than expected 
for quasi-static equilibrium and live for much longer than a cooling time. The loops are 
therefore neither steadily heated nor cooling as monolithic structures. It has been shown 
that the observed densities and timescales can be explained by bundles of nanoflare he ated 
strands, as long as nanoflares do not all occur at the same time (see iKlimchukl 120061 and 
references cited therein). Neighboring strands will therefore have different temperatures, 
and loops are predicted to have multi-thermal cross sections. In particular, emission should 
be produced at temperatures hotter than 3 MK. Hot loops are sometimes observed at the 
locations of warm loops, but not always. Example 5 suggests that nonthermal electron 
beams are a possible explanation for the lack of hot emission. As we have discussed, beams 
can produce excess densities through evaporation without the need for high temperatures. 
We have just begun to explore this possibility. For now, we note that the coronal DEM(T) 
curve of example 5 ( Figure [8]) bears a close resemblance to the observed curves reported by 
Schmelz et al.l (120011 ) for a loop seen above the limb. 



In conclusion, EBTEL is a powerful new tool that can be applied to a variety of prob- 
lems in which large numbers of evolving strands must be computed. For example, it is now 
feasible to construct multiple models of nanoflare-heated active regions or entire stars and 
therefore to examine a wide array of nanoflare parameters (magnitude, lifetime, occurrence 
rate, dependence on field strength and strand length, etc.). By determining which parame- 
ters best reproduce the observations, we can place import ant constraints on the heating and 
thereby gain insight into the physical mechanism (e.g., iMandrini. Demoulin. fc Klimchuk 
2000l ; ISchrijver et al.l 12004 ; IWarren fc Winebargerl 120061 ) . EBTEL is currently being use d 



to study the emission characteristics of coro nal arcades (IPatsourakos fc K 



imchukl l2007bl ). 
2007J) . and to 



to explain the light curves of solar flares (IRaftery. Gallagher. &: Milliganl . 
model coronal loops as self-organized critical systems (jLopez Fuentes &: Klimchuk 2005 
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Klimchuk. Lopez Fuentes. fc DeVordl2006l ). 



Interested users are invited to contact us for a copy of our IDL code. 
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providing comparison ID simulations. The authors benefited from participation on the 
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understanding coronal heating (Team Parenti). 



A. Transition Region DEM{T): Alternate Derivation 

An alternate approach to deriving the differential emission measure distribution of the 
transition region is to consider three limiting cases — strong evaporation, strong condensation, 
and static equilibrium — and to combine the results into a single form with smooth transitions. 



A.l. Strong Evaporation 

During strong evaporation, the heat flux from the corona far exceeds the radiative losses 
from the transition region, |F | 3> TZtr, and the energy equation reduces to an approximate 
balance between thermal conduction heating and enthalpy cooling: 

We here use the classical form for the heat flux because saturation is not expected with the 
relatively low temperatures and high densities of the transition region. Integrating equation 
( 1Aip . we obtain 

f - ^ < A2 > 

where we have used the ideal gas law and equation for the constant electron flux. Sub- 
stituting into equation (1281) and noting that pressure is approximately constant throughout 
the transition region and corona, we have the final expression 

DEM e JT) « — ^— (A3) 
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A. 2. Strong Condensation 

During strong condensation, the heat flux from the corona is much less than the radiative 
losses from the transition region, \Fq\ <C IZtr, and the energy balance is then between enthalpy 
heating radiation cooling: 

n 2 A(T) » -~{Pv). (A4) 
2 as 

The constant electron flux allows us to write 

= 2*J„f. (A5) 

SO 

dT n 2 A(T) 

— ss — (A6) 

ds 5kJ 1 ; 

and 

DEM con (T) « -^y- (A7) 



A. 3. Static Equilibrium 

Lastly, in static equilibrium, the heat flux from the corona very nearly balances the 
radiative losses from the transition region, |F | ~ 1Z tr - The inequality is broken only by a 
possible source of direct plasma heating, which is likely to be very small in comparison to 
thermal conduction heating. Nonthermal electrons are a possible exception. Barring this 
possibility, 



2 T 7 / 2 



^o— , (A9) 



where 



H * = MISi (A10) 



is the temperature scale height. This gives 

dT ( 7 \ 1,2 nk(Tfl 2 



ds \2k J T 3 / 4 

and 



(All) 



M Mr).(g) 1/2 M(r ^ rl/ , (Ai2) 
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We can combine these three limiting cases into a single expression that applies at all 
times: 

DEM{T) = (f DEM £V - f° n ^ DEM se + 1Z tr DEM con ) ( F - + K tr 

V -^0 + 'Mr / \ -TO + 'Mr 

(A13) 

This expression reduces to the desired forms in the relevant limits. We have confirmed 
that the temperature dependencies in equations (IA3I) . flA7j) . and flA12j) . are present in the 
differential emission measure distributions from ID simulations. 
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Table 1: ID Equilibrium Parameters (Short Strand) 



T a (MK) 




f/T a 


T /T a 


0.80 


2.5 


0.88 


0.57 


1.83 


1.8 


0.89 


0.61 


3.77 


6.7 


0.87 


0.46 


4.60 


9.5 


0.87 


0.40 


7.08 


17.1 


0.86 


0.28 


10.40 


20.7 


0.86 


0.22 
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Table 2: ID Equilibrium Parameters (Long Strand) 



T a (MK) 




f/T a 


T /T a 


1.00 


4.7 


0.90 


0.59 


1.94 


2.9 


0.90 


0.63 


3.95 


4.3 


0.90 


0.61 
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Example 1 
51 ' ' ^ 




oi 



O 1000 2000 3000 4000 5000 

Time (s) 



^ 6 IT 

I 




1000 2000 3000 4000 5000 

Time (s) 





0.30 : 


CM 
I 


0.25 r 


E 


0.20 [ 


O 

c 


0.15 [ 



"D 



0.10 r 

0.05 'r. 

0.00 K' 





1000 2000 3000 
Time (s) 



4000 5000 



Fig. 1. — Evolution the coronal- averaged temperature, electron density, and pressure for 
a loop strand heated impulsively by a 500 s nanoflare (example 1). Solid curves are for 
the EBTEL simulation, and dashed curves are for the ID simulation. Classical heat flux is 
assumed at all times. 
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Fig. 2. — Differential emission measure distribution for the whole strand (unit cross section) 
averaged over the first 10 4 s of the 500 s nanoflare simulation (example 1). Solid curve is for 
the EBTEL simulation, and dashed curve is for the ID simulation. 
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Fig. 3. — DEMi n (T) = T x DEM(T) corresponding to the differential emission measure 
distributions in Fig. [2j 



22 



Example 1 




Fig. 4. — Coronal (dashed) and transition region (dot-dashed) contributions to the total 
differential emission measure distribution (solid) from the EBTEL simulation of the 500 s 
nanoflare (example 1), and coronal contribution from the ID simulation (dotted). 
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Fig. 5. — Evolution the coronal- averaged temperature, electron density, and pressure for a 
loop heated impulsively by a 50 s nanoflare (example 2). Solid curves are for the EBTEL 
simulation, dashed curves are for the ID simulation, and dotted curves are for the Cargill 
simulation. Classical heat flux is assumed at all times. 
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Fig. 6. — Evolution the coronal- averaged temperature and electron density for the time- 
dependent coronal heating rate shown at the bottom (example 3). Solid curves are for 
the EBTEL simulation and dashed curves are for the ID simulation. Classical heat flux is 
assumed at all times. 
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Fig. 7. — Total (solid), coronal (dashed), and transition region (dot-dashed) differential 
emission measure distributions for an EBTEL simulation of a nanoflare that is ten times 
larger than that of example 1. 
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Fig. 8. — Coronal differential emission measure distribution for example 1 (dashed) and 
for a corresponding simulation in which half of the nanoflare energy takes the form of a 
nonthermal electron beam (solid). 



